Etude prévisionnelle des ventes

Abdelghafour Ait bennassar, Salma Mounji

Table of Contents

Importation des outils nécessaires

Analyse exploratoire des données

Importation de données

Description de données

Initialement, on disposait des données sous forme d’un classeur Excel contenant cinq feuilles où chaque feuille contient les ventes mensuelles des produits dans un an, pour obtenir les données sous la forme requise pour l’étude, on a joint les données des cinq années pour obtenir un jeu de données dans lequel les lignes représentent les dates (63 mois) et les colonnes représentent les produits (70 produits).

Le corrélogramme : visualisation de la matrice de corrélation

On observe une forte correlation entre certains groupes de produits, ce qui suggère qu'on peut les classifier selon leurs comportements de ventes.

clustering : Regroupement des produits

Normalisation des variables

L'un des problèmes qu'on doit traiter avant de commencer le partitionnement (clustering) est l'échelle de la série. Sans normaliser les données, les séries qui se ressemblent seront vues si différentes les unes des autres et affecteront la précision du processus de clustering.
Nous pouvons voir l'effet de la normalisation dans les images suivantes.

Réduction de la dimension : Analyse en composantes principales (ACP)

Afin de regrouper nos séries avec des k-moyennes, les indices temporels des séries chronologiques seront considérés comme des caractéristiques différentes et seront les dimensions des points de données(les séries).
Puisque la longueur des série temporelles en pratique est souvent importante, un autre problème auquel on doit faire face est la malédiction de la dimensionnalité,ce terme a été inventé pour la première fois par Richard E. Bellman lors de l'examen des problèmes de programmation dynamique. Cela signifie essentiellement que lorsque la dimensionnalité des données augmente, la distance entre les points de données augmente également. Ainsi, ce changement de mesure de la distance affecte gravement les algorithmes basés sur la distance.
Pour résoudre ce problème, on va effectuer une analyse en composantes principales ACP avant le partitionnement.

Désormais avec moins de dimensions qu'avant, nous pouvons voir comment nos séries se répartissent en 2 dimensions.

Clustering : K-means

Le partitionnement en k-moyennes (ou k-means en anglais) est une méthode de partitionnement de données et un problème d'optimisation combinatoire. Étant donnés des points et un entier k, le problème est de diviser les points en k groupes, souvent appelés clusters, de façon à minimiser une certaine fonction. On considère la distance d'un point à la moyenne des points de son cluster (centroid) ; la fonction à minimiser est la somme des carrés de ces distances,la distance euclidienne est souvent utilisée.

Initialement on va éstimer le nombre de clusters requis par la racine du nombre des points (séries), on obtient 9 clusters.

et on va appliquer la méthode de Kmeans sur le résultat de l'ACP.

Après l'entrainement du modèle, on a présenté les résultats. Pour chaque cluster, on a tracé les séries en gris, et afin de visualiser le mouvement et le comportement des séries appartenant au même cluster, on a pris la moyenne du cluster et tracé cette série moyenne en rouge.

Nous pouvons voir la distribution des séries chronologiques en clusters dans le graphique suivant.
les séries semblent uniformement distribuées aux clusters.

Pour bénéficier de partitionnement qu'on a fait, il est nécessaire de lister pour chaque cluster, les séries qui lui appartiennent.

L’un des avantages que nous offre le regroupement effectué est l’exploitation de l’information disponible non seulement sur les ventes du produit qu’on souhaite prévoir, mais également l’information sur les produits ayant un comportement de ventes similaire (produits appartenant au même groupe).

On a selectionné le produit le plus proche du centroide de chaque cluster, autrement dit le plus proche aux autres produits au même groupe.

Pour visualiser ce résultat, on a tracé l'ensemble des points ainsi que les centres des groupes (étoiles), le point le plus proche de chaque centre est présenté en gris.

pour re vérifier ces résultats, on peut voir que les séries choisis sont centrées pour chaque groupe.

Analyse de données

On a extrait les produits représentatifs des groupes et on a construit un nouveau jeu de données de séries temporelles sur lequel on va travailler par la suite.

Statistiques descriptives :

Dans cette section, nous allons trouver quelques mesures de tendance et de dispersion comme la moyenne, médiane, la variance et l’écart type pour mieux percevoir les caractéristiques des médicaments.

D'après le tableau ci-dessus et le graphique, il est clair que ces médicaments ont des natures diverses et nous pouvons étendre cette remarque à tous les médicaments. Leurs moyennes et leur somme de ventes sont très différentes. Certains d'entre eux ont été vendus fortement comme le médicament 30, et d’autres sont moins vendus comme les médicaments 47 et 64. Certains médicaments ont une variance élevée comme les médicaments 30 et 22. Par conséquent, nous pouvons conclure que la majorité des médicaments ont des caractéristiques différentes en plus de la dissemblance dans le comportement de vente.

Visualisation de la tendance et saisonnalité

On peut observer une saisonnalité dans certains produits notamment P_64 et p_47, avec une variété de tendances, certaines gardent la même monotonie au cours du temps tel que celles des P_64 et d’autre produits possèdent une tendance qui alterne entre haussière et baissière.
L’existence de hausses soudaine et significative dans le produit P_22 pourra biaiser la prévision.

Pour mieux visualiser la saisonnalité et tendance des ventes des médicaments, on a tracé pour chaque produit, une moyenne mobile 4-mois (pour détecter la saisonnalité) ainsi qu’une moyenne mobile 12-mois (pour capter la tendance)

Les graphiques ci-dessous représentent la décomposition des séries temporelles

L'analyse des portions de chaque composant de la série chronologique est particulièrement utile pour déterminer l'absorption des résidus dans les données, sur la base des données décomposées. Le volume de cette absorption implique la prévisibilité de la série chronologique - plus les résidus sont élevés, moins la prévisibilité.

Le tableau ci-dessous présente la moyenne des observations et la moyenne des résidus pour chaque produit, ainsi que la proportion de l’erreur par rapport à la série chronologique observée, pour l’ensemble des médicaments l’erreur ne dépasse pas 7% des ventes observées, chose qui indique qu’une quantité importante d’information pourra être expliquée.

Modélisation statistique univarié : ARIMA

Analyse de la stationnarité

Analyse d’autocorrélation

L'analyse d'autocorrélation illustre le potentiel de prédiction des données de séries chronologiques. Les graphiques d'autocorrélation résument graphiquement la force d'une relation avec une observation dans une série chronologique avec des observations à des pas de temps antérieurs. Le coefficient de Pearson est utilisé pour mesurer l'autocorrélation. Ainsi, l'analyse suivante n'est pertinente que pour les données avec une distribution gaussienne normale.

Un tracé de l'autocorrélation d'une série chronologique par décalage est appelé la fonction d'autocorrélation (ACF). Ce graphique est parfois appelé corrélogramme ou graphique d'autocorrélation. Le graphique montre la valeur de décalage le long de l'axe des x et la corrélation sur l'axe des y entre -1 et 1. Les intervalles de confiance sont dessinés sous forme de cône. Par défaut, il est défini sur un intervalle de confiance de 95%, ce qui suggère que les valeurs de corrélation en dehors de ce code sont très probablement une corrélation.

La figure ci-dessous présente la fonction d’autocorrélation (ACF) pour chacun des médicaments avec un intervalle de confiance de 95%.

D'après les correlogrames des séries temporelles ci dessus, on peut constater que la fonction d'autocorrelation se dégrade ou diminue lentement, ce qui indique la non stationnarité des séries.
De plus, dans les deux séries des produits p_47 et p_64 l'ACF montre une oscillation et des pics qui se reproduisent avec des décalages de 12 mois indicant un effet de saisonnalité.

Tests de stationnarité

Afin d'estimer le terme de différenciation saisonnière D, on a effectué deux tests de saisonnalité pour différents niveaux de D pour estimer le nombre de différences saisonnières nécessaires pour rendre la série chronologique stationnaire et on a sélectionné la valeur maximale de D pour laquelle la série chronologique est jugée saisonnièrement stationnaire par le test statistique.

De même, pour estimer le terme de différenciation d, on a effectué trois test statistiques de racine unitaire : Augmented Dickey Fuller (ADF), Phillips–Perron (PP) et Kwiatkowski–Phillips–Schmidt–Shin (KPSS) pour différents niveaux de d afin d'estimer le nombre de différences nécessaires pour rendre chaque série chronologique stationnaire, et on a sélectionné la valeur maximale de d pour laquelle la série temporelle est jugée stationnaire par le test statistique.

Le tableau ci dessous synthétise les résultats, chaque test de saisonnalité indique la présence ou absence d'un effet de saisonnalié, et de même pour les tests de racine unitaire, ainsi que les ordres de differenciation nécaissaires pour stationnarier la série.

Stationnarisation des séries

On a procedé par la suite à la stationnarisation des séries en appliquant la differenciation adéquate, et on a reéffectué les mêmes tests pour confirmer qu'on a atteint la stationnarité des séries.

Les tests finals confirment qu’on a bien atteint la stationnarité des séries.

Prévisions par ARIMA

selection des meilleurs modèles

Les graphiques ci-dessous représentent pour chaque produit la fonction d’autocorrelation (ACF) ainsi que la fonction d’autocorrélation partielle (PACF) de la série stationnarisée.

A l’addition de la stationnarité des séries qui peut être confirmé à partir de ces corrélogrammes, ces derniers nous permettent aussi de déterminer l'ordre d'ARIMA, mais la détection visuelle ne permet pas toujours d'atteindre la performance maximale, c'est pourquoi on a utilisé la fonction Auto_Arima qui permet d'ajuste le meilleur modèle ARIMA à une série temporelle univariée selon un critère d'information fourni (soit AIC, AICc, BIC...). La fonction effectue une recherche (soit pas à pas, soit parallélisée) sur des ordres de modèle et saisonniers possibles dans les contraintes fournies, et sélectionne les paramètres qui minimisent la métrique donnée.# L'approche parallèle est une méthode de brute force naïve, il s'agit d'une recherche de grille (grid search) sur diverses combinaisons d'hyperparamètres. Cela prendra généralement plus de temps pour plusieurs raisons. Tout d'abord, il n'y a pas de procédure intelligente sur la façon dont les commandes de modèles sont testées ; ils sont tous testés (pas de court-circuit), ce qui peut prendre un certain temps.

L'approche par étapes (stepwise) suit la stratégie définie par Hyndman et Khandakar dans leur article de 2008, « Automatic Time Series Forecasting : The Forecast Package for R ».

Étape 1 : Essayez quatre modèles possibles pour commencer :

ARIMA(2,d,2) si m = 1 et ARIMA(2,d,2)(1,D,1) si m > 1
ARIMA(0,d,0) si m = 1 et ARIMA(0,d,0)(0,D,0) si m > 1
ARIMA(1,d,0) si m = 1 et ARIMA(1,d,0)(1,D,0) si m > 1
ARIMA(0,d,1) si m = 1 et ARIMA(0,d,1)(0,D,1) si m > 1
Le modèle avec le plus petit AIC (ou BIC, ou AICc, etc., selon les critères de minimisation) est sélectionné. C'est le "meilleur" modèle actuel.

Étape 2 : Considérez un certain nombre d'autres modèles :

Où l'un des p, q, P et Q peut varier de ±1 par rapport au meilleur modèle actuel
Où p et q varient tous les deux de ±1 par rapport au meilleur modèle actuel
Où P et Q varient tous les deux de ±1 par rapport au meilleur modèle actuel
Chaque fois qu'un modèle avec un critère d'information inférieur est trouvé, il devient le nouveau meilleur modèle actuel, et la procédure est répétée jusqu'à ce qu'il ne puisse pas trouver un modèle proche du meilleur modèle actuel avec un critère d'information inférieur ou que le processus dépasse l'un des seuils d'exécution défini via arima.StepwiseContext.

Les meilleurs modèles sélectionnés sont :

Prévisions statiques et dynamiques

On a effectué la prévision de l'ensemble des séries par les meilleurs modèles ARIMA séléctionnés selon deux approches.

La méthode de prévision statique calcule une séquence de prévisions à un pas d'avance, en utilisant les valeurs réelles plutôt que prévues pour les variables dépendantes décalées, si elles sont disponibles. En effet, pour effectuer la prévision d'un nouveau mois, le modèle est mis a jour en incluant la valeur réelle du mois précédant, cette méthode est plus performante et donne des prévisions moins biaisées mais ne peut etre utilisée que pour prévoir un seul mois à l'avance.

La méthode de prévision dynamique calcule des prévisions dynamiques en plusieurs étapes à partir de la première période de l'échantillon de prévision, les valeurs précedemment prévues pour les variables dépendantes décalées sont utilisées pour former la prévision de la valeur actuelle. Cette méthode est utile lorsque l'entreprise s'interesse à effectuer une prévision à long terme, c'est dire prévoir les ventes de plusieurs mois à l'avance.

Les graphiques ci-dessous représentent les prévisions effectuées par les deux méthodes ainsi que les intervalles de confiance (95%). On peut constater que les séries des deux produits P_22 et P_30 correspondent à des marches aléatoires, la meilleure prévision dans ce cas est la valeur précendente. Les intervalles de prévision associés à la méthode dynamique s'allongent à mesure que l'horizon de prévision s'allonge, plus nous prévoyons à l'avance, plus l'incertitude est associée à la prévision, et donc plus les intervalles de prévision sont larges. Dans les prévisions par la méthode statique, la variation de la distribution des prévisions est presque la même que la variation des résidus.

Diagnostic des résidus

Pour la validation des modèles obtenus, on a effectué un diagnostic des résidus. le tableau ci-dessous synthétise les résultats des tests de bruit blanc, homoscedasticité et normalité pour le résidus de chacune des séries

Modélisation univarié par l’approche d’apprentissage automatique

Features engineering

Pour la démonstration de la demarche suivie dans l'apprentissage supervisé, on a selectionné aléatoirement un produit (P_64).

On a commencé par implémenter le concept de "feature engineering" ou on a généré les variables nécessaires pour entrainer les modèles d'apprentissage.

Pour résoudre ce problème, on a effectué une transformation par la fonction cosinus après avoir normalisé les variables entre 0 et 2π, ce qui correspond à un cycle de cosinus, cette solution ne régle pas le problème complétement car deux valeurs différentes peuvent avoir la même image par la fonction cosinus. La meilleure façon de résoudre ce nouveau problème était d'ajouter une autre information cyclique pour distinguer deux temps avec des valeurs de cosinus identiques, il s'agit de la fonction sinus. Nous pourrions le considérer comme un système de coordonnées à deux axes.

Comme indiqué précédemment, on a utilisé les variables de la série décalées pour beneficier de l'autocorrélation existante dans la série. Pour prendre en compte l'effet de la saisonnalité, on a considéré 12 retards.

On a généré les moyennes mobiles et ecart-types mobiles (trimestrielle, semestrielle et annuelle) des valeurs décalées afin de détecter la tendance globale des données.

Modèles de régression

Régression linéaire multiple

On a ajusté un modèle de regression linéaire multiple sur l'ensemble des variables.

On peut observer que la courbe de régression est collée à la série temporelle du produit, ce qui indique qu’il s’agit d’un modèle complexe et flexible qui surapprend les données d’entrainement (un biais faible) mais conduit à des erreurs plus importantes sur les nouvelles données (une variance élevée). La figure de l’importance des variables confirme le surapprentissage. En effet, le modèle concentre toutes ses forces sur les variables du temps et néglige l’information inclue dans les autres variables. Lorsque les variables indépendantes sont corrélées, cela indique que les changements dans une variable sont associés à des changements dans une autre variable. Plus la corrélation est forte, plus il est difficile de changer une variable sans en changer une autre. Il devient difficile pour le modèle d'estimer la relation entre chaque variable indépendante et la variable dépendante indépendamment parce que les variables indépendantes ont tendance à changer à l'unisson. Dans notre cas, plusieurs variables sont fortement corrélées comme l’indique le corrélogramme suivant :

D'autre part, puisqu’on a différentes échelles dans les variables, celles-ci doivent être transformées en une même échelle pour explorer l'importance des caractéristiques et être utilisés par la suite dans la régularisation.

modèles de régularisation : Ridge et Lasso

Pour une meilleure optimisation des variables explicatives, on a appliqué la régularisation. Deux des modèles de régression avec régularisation les plus populaires sont les régressions Ridge et Lasso. Ils ajoutent tous deux des contraintes supplémentaires à notre fonction de perte.

Dans le cas de la régression Ridge, ces contraintes sont la somme des carrés des coefficients multipliée par le coefficient de régularisation. Plus le coefficient d'une caractéristique est grand, plus notre perte sera importante. Par conséquent, nous essaierons d'optimiser le modèle tout en gardant les coefficients assez bas.

À la suite de cette régularisation L2, nous aurons un biais plus élevé et une variance plus faible, donc le modèle se généralisera mieux.

Le deuxième modèle de régression, la régression Lasso, ajoute à la fonction de perte, non pas des carrés, mais des valeurs absolues des coefficients. Par conséquent, au cours du processus d'optimisation, les coefficients des caractéristiques sans importance peuvent devenir des zéros, ce qui permet une sélection automatisée des caractéristiques. Ce type de régularisation est appelé L1. Les prévisions par les deux régularisations ainsi que le graphique de l’importance des variables sont présentées dans les figures suivantes :

On peut remarquer que la qualité d’ajustement s’est améliorée par la régression Ridge, l’erreur sur le train set est la même que l’erreur sur le test set, donc le modèle ne surapprend pas les données d’entrainement et le problème de surapprentissage est réglé. On peut également remarquer que les intervalles de prédiction sont moins larges que dans la régression linéaire. Nous pouvons clairement voir que certains coefficients se rapprochent de plus en plus de zéro (bien qu'ils ne l'atteignent jamais réellement) à mesure que leur importance dans le modèle diminue.

De même que la régression Ridge, la régression Lasso permet de régler le problème de surapprentissage. La régression au lasso s'est avérée plus conservatrice, le modèle elle a supprimé plusieurs variables.

La figure ci-dessous représente les prévisions pour chaque produit par le modèle qui minimise l'erreur (MAPE) parmis les deux modèles de régularisation.

Gradient Boosting : XGBoost

On a effectué les prévisions par un autre modèle d’apprentissage, à savoir XGBoost. Ce modèle est très performant sur beaucoup de problèmes de prédiction, évaluons sa performance sur notre problème de prévision des séries temporelles.

L’implémentation XGBoost du grandient Boosting propose plusieurs hyperparamètres tel que :
-Max_depth : Profondeur maximale d'un arbre. L'augmentation de cette valeur rendra le modèle plus complexe et plus susceptible de surapprentissage.
-colsample_bytree est le sous-échantillon des colonnes considéré lors de la construction de chaque arbre. Le sous-échantillonnage a lieu une fois pour chaque arbre construit.
-max_leaves :Nombre maximum de nœuds à ajouter.
On considère encore une fois le produit P_64, on a essayé plusieurs combinaisons des hyperparamètres jusqu’à trouver la combinaison optimale qui minimise l’erreur (MAPE), on a entrainé le modèle en utilisant cette dernière, les prévisions données sont présentées dans le graphique ci-dessous.

On a réeffectué le tunning pour chacun des produits, choisi la combinaison des hyperparamètres qui minimise l’erreur et ajusté un modèle adapté. Le tableau ci-dessous présente l’erreur des modèles entrainés et testé sur le test set pour chaque produit.

Finalement, pour chaque produit on a retenu le modèle qui minimise l’erreur sur le test set parmi les modèles d’apprentissage supervisé testés : RLM, Ridge, Lasso, XGBoost. Les prévisions données par le modèle qui minimise l’erreur pour chaque produit ainsi que son erreur sur le test set sont présentées dans le graphique suivant.

Comparaison des modèles univariés

Modélisation multivariée : Vector AutoRegressif Moving Average (VARMA)

Dans cette section, on a exploité le clustering effectué auparavant en s’intéressant à l’ensemble des produits d’un même groupe (Cluster).
Ayant des comportements de ventes similaires, ces produits peuvent s'influencer mutuellement et par la suite l’information supplémentaire sur les évolutions des autres produits peut améliorer la qualité du modèle.
On dispose de 9 clusters, Pour la démonstration on considère le cluster qui comporte les produits :P_3, P_21, P_23, P_24, P_25, P_27. Visuellement, on peut constater que les produits ont des évolutions similaires avec des effets de saisonnalité semblables.

Le corrélogramme suivant confirme la forte corrélation entre les produits, plus particulièrement les produits P_21, P_23 et P_24.

D’une manière analogue à ARIMA, le modèle VARMA nécessite la stationnarité des séries temporelle. Pour vérifier la stationnarité, on a effectué les mêmes tests de racine unitaire et de saisonnalité, le tableau suivant synthétise les résultats.

Après les différenciations nécessaires, on obtient les séries suivantes :

Pour trouver le meilleur modèle VARMA, on a utilisé le module AutoTS, qui fonctionne actuellement de la manière suivante :

Le processus commence par le remodelage et le pré-traitement de données selon le besoin.
Un premier fractionnement train/test est généré où le test est de longueur 'Forecast_length'.
Le modèle initial est une combinaison d'apprentissage par transfert et de modèles générés aléatoirement. Ceci est testé sur le train/test initial.
Les modèles se composent d'une étape de pré-transformation (options de remplissage, options de suppression des valeurs aberrantes, etc.)
Les meilleurs modèles (sélectionnés par une combinaison de métriques) sont recombinés avec des mutations aléatoires pour n_générations.
Un pourcentage des meilleurs modèles issus de ce processus est soumis à une validation croisée, où ils sont réévalués sur de nouvelles divisions train/test.
S'il est utilisé, l'assemblage horizontal utilise les données de validation pour choisir le meilleur modèle pour chaque série.
Le meilleur modèle ou ensemble en cours de validation est sélectionné comme best_model et utilisé dans la méthode .predict() pour générer des prévisions.

Les prévisions par le modèle selectionné sur les séries donne les résultats suivants :